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Abstract 

A model describing the dynamics of crystal exciton states on a 
hexagonal lattice has been studied in [1]. The hamiltonian includes 
terms describing the movement of the exciton throughout the crystal. 
Here we generalize the hamiltonian by including first order couplings 
of such terms to the lattice phonon fields. A discussion of the ground 
state configuration of the exciton field as a function of the system 
couplings is given. Finally we analyse the stability of certain excited 
states under perturbations. 

1 Introduction 

The promotion into the conduction band of an electron in a solid results in 
the creation of a hole in the valence band, which can be viewed as carrying a 
positive charge. The electron and hole therefore feel a Coulomb attraction. 
A quasiparticle may result, known as an exciton, with an energy lower than 
the bottom of the conduction band. The conduction electron acts to partially 
screen the interaction between the surrounding nuclei and this results in a 
modification to the lattice potential. The neighbouring sites are drawn in 
and this produces a 'self-trapping' of the exciton. 

The importance of this property of the exciton was first realised by Davy- 
dov and Kislukha in 1973 [2]. Working in one dimension, they proposed that 
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the exciton state could explain the dispersionless transport of energy in pro- 
tein chains within the human body - a problem that had been important 
enough to earn the title: 'crisis in bioenergetics' [3]. The quasiparticle moves 
along the chain, accompanied by its associated lattice deformation. Due to 
its solitary wave-like properties, the bound state was termed a soliton. The 
existence of this state was found to exhibit a crucial dependence upon the 
electron-phonon coupling strength. 

The system used in modelling the one-dimensional chain was more re- 
cently adapted to describe the two-dimensional case. In the simplest exam- 
ple, i.e. the square lattice, an analysis of the stationary case showed that 
in the continuum limit the equations governing the evolution of the system 
could be reduced to a single Nonlinear Schroedinger Equation with an ex- 
tra term resulting from the discreteness of the lattice [4-6]. This equation, 
without the extra term, is well known to possess localized solutions in one 
dimension [7]. In two dimensions however, the soliton is stable only in the 
presence of the extra term. 

The equations were then adapted to the case of a hexagonal lattice 
[1]. Again sohton solutions were found to exist. Furthermore, the discrete 
equations were found, in the stationary limit, to obey a Discrete Nonlinear 
Schroedinger Equation. A similar equation was the starting point for inves- 
tigations on the hexagonal and triangular lattices [8] . The fact that no such 
discrete solution could be found in the case of a square lattice suggests that 
solitons could have an important role in the dynamics of hexagonal crystal 
lattices, and therefore carbon nanotubes. 

Nanotubes can be thought of geometrically as rolling up a two dimensional 
hexagonal graphite sheet. The angle at which the sheet is rolled, known 
as the chiral angle [9], determines whether the carbon nanotube is semi- 
conducting or metallic. Because of the band properties of semiconductors, 
it is these types of tubes in which cxcitons are most important. Due to the 
nanoscale structure of carbon nanotubes, even a single exciton state could 
have a considerable effect on the structure's electrical properties. 

The system Hamiltonian includes terms that describe the movement of 
an exciton to a neighbouring lattice site, known as the resonance interaction. 
Studies so far have worked in the zeroth approximation excluding terms cou- 
pling the resonance interaction to the lattice. Here we present the results 
of our investigations taking into account first order couplings of the reso- 
nance interaction to the lattice phonons: terms expected to be important in 
describing carbon naotubes. 
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2 The Hamiltonian and Equations of Motion 



Our system is described by a Hamiltonian containing three parts, the first 
part corresponding to the original electron Hamiltonian discussed in [1] (in- 
cluding interaction terms with the lattice), the second to the lattice potential 
and the third to the new interaction terms introduced in this paper. 



H = He + Hph + Hint (1) 

There are two distinct types of site within the hexagonal lattice, which 
must be summed over separately in the Hamiltonian. Denoting the distinct 
types imaginatively as type 1 and type 2 (see figure 1), the Hamiltonian is: 




Figure 1: The two inequivalent types of site in the hexagonal lattice. 
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The bosonic creation and annihilation operators for the exciton at site 
are given by Ajj and Aj respectively. The operators iiij and Vij give 
the displacement of molecule in the x and y directions, pij and qij give 
their respective conjugated momenta. The zeroth order coupling in the res- 
onance interaction of the exciton is denoted j, the new term describing first 
order couplings of the resonance interaction to the lattice has a strength /. 
The coupling constant k describes the harmonic lattice potential, which, in 
the limit of small amplitude vibrations, is assumed to provide a good approxi- 
mation to the true crystal potential. Terms of first order in the displacements 
in Hph are zero due to the associated coupling vanishing in equilibrium. The 
exciton part of the original Hamiltonian Hg contains the coupling c, which 
describes the displacement of the lattice in the presence of an exciton. Physi- 
cally the presence of an exciton at the site changes the potential between 
neighbouring molecules by modifying the interatomic Coulomb force, since 
an electron has been promoted into a more weakly bound conduction band 
state. The modification to the Hamiltonian involving the addition of terms 
proportional to / is given by Hint- 
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The analysis proceeds in the usual way [4-6] by allowing the Hamiltonian 
to operate on the approximate product wavefunction: 



(5) 



where is an excited state of the lattice and \^{t)) is a coherent phonon 

state. The states A1j|0) are complete and orthonormal. The ipij are the 
complex coefficients of: 



(6) 



The displacement and momentum operators act upon |$) to give: 

=p,j; = qi,j 
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The semiclassical Hamiltonian can then be written in terms of the complex 
amplitudes ipijit) by calculating: 
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which tranforms the hamiltonian to: 
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where the phonon energy W is given by: 
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The Hamiltonian is thus completely written in terms of the complex field 
ij^ij and the actual molecular displacement fields Uij and Vij. The full equa- 
tions of motion for the Hamiltonian (9) are given in Appendix Al. Below 
we give the equations of motion corresponding to the case when the lattice 
site is of type 1, after performing a set of rescalings to reduce the number of 
parameters. 
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where the following rescalings have been performed: 
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We have attempted to find solutions to these equations in the stationary 
limit, both for the discrete case and in the continuum limit. Unfortunately 
no simple solution was found. In [1] it was found that one could express the 
discrete four-point laplacian of the relevant combination of the displacement 
fields in terms of a seven-point laplacian in the exciton field. The hexago- 
nal lattice admitted a simple solution to this equation and this allowed the 
expression of the system equations of motion in terms of a single discrete 
nonlinear Schroedinger equation with an extra stabilizing term. The addi- 
tion of the terms in / makes such a task much more difficult, and as yet no 
such solution has been found. 



3 Numerical Simulation of the Equations 

Simulation of the equations (11-13) involved the use of a fourth order Runge- 
Kutta method. In order to reduce the simulation time a compression of the 
lattice was performed, as discussed in [1]. The Hamiltonian and equations 
of motion that were directly simulated are presented in Appendix Al. In 
order to represent carbon nanotubes to a first approximation, continuous 
boundary conditions were imposed upon the lattice, in both directions. The 
longitudinal lattice dimension was chosen to be significantly larger than the 
transverse direction. The lattice size was set to 20 x 8 sites. Although this 
was a relatively small lattice size, it enabled a detailed invcstiagtion of the 
system states. A future aim should be to investigate the states on a larger 
lattice. 
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3.1 Finding the Minimum Energy 



3.1.1 Absorption of Kinetic Energy 

In order to find and then investigate the static configuration through simu- 
lation it was necessary to absorb energy from the equations. One method of 
achieving this involved the addition of the following damping terms to the 
respective right-hand sides of (4) and (5): 



and -u^, (15) 

where ly determines the strength of the damping. These terms have the effect 
of gradually absorbing the kinetic energy of the molecules vibrating about 
their equilibrium positions. Figure 2 shows the descent of the system energy 
for three different values of the damping parameter u. The parameters used 
for simulation were A; = l,j = l,/ = and c = 23.4. The initial conditions 
involved field configurations corresponding to a polaronic like state, i.e. a 
state localised mainly around a single lattice site. It was known that this 
state was the minimum energy state for the values of the couplings stated 
above. The damping was then used to reduce the system energy towards 
the static configuration. In this particular case it was found that u ^ 1.0 
produced the quickest descent to the minimum energy, however varying the 
couplings has an effect upon the optimum damping coefficient. Similar effects 
were observed for other values of / and c. 



3.1.2 Absorbtion of Energy from the Exciton Field 

The method described in section 3.1.1 was effective in reducing the system 
energy. However, it was found that the speed of the process could be greatly 
increased by introducing a direct absorption of energy from the exciton field. 
The equation of motion for t/^ (11) contains only a first order derivative with 
respect to time, and so the introduction of a damping term similar to that 
in (15) would destroy the system dynamics. An approximate method was 
therefore introduced so that the damping term could be included. 

For a quantum state with a definite energy E, the time dependence of the 
wave function can be taken to be that of a plane wave. The eigenfunctions 
(t){x) are then the stationary states of the system: 

ipi,j{x, t) = 4)ij{x)exp(^-^^^ (16) 
The variation of i/j over an infinitessimal time step t is then: 



8 



lime 



■4 



t> 1 i i ^ ^ 1 1 1 1 


1 1 a s J 1 1 1. 1 1 


L < 


It 

\\ ^ 


\\ 




i g 









Figure 2: Variation of the system energy due to the addition to the equations 
of motion of the terms in (15). Three different damping parameters v are 
shown. 



In the presence of damping, the stationary states t) have some small 
time dependence. The value of the energy E was known throughout simula- 
tion. The variation of i/j given by (16) could therefore be taken into account. 
It was then possible to introduce a damping term that would give the varia- 
tion dt(f)ij{x, t) in (17) without affecting the variation of the stationary states 
according to (16). 

The LHS of (17) is specified by the equation of motion (11). Using (16) 
to re-write ip in terms of (j), and cancelling the exponential functions, yields 
the equation: 



ih^^-E(l),, + F{cP,,,), 
9 



(18) 



where Fi^cpij) is the RHS of equation (11) with replacing ijj. The right 
hand side of this equation was calculated within the simulation. A damping 
term could then be introduced without destroying the system dynamics: 



{a + ^h)^ = -E<t>^, + F{<t>^,) (19) 

where the magnitude of a specifies the strength of the absorption. 

The implementation of this method resulted in an enormous reduction 
in the time taken to reach the stationary configuration. It was found that 
a combination of the two methods for extracting the energy (15) and (19) 
produced the quickest descent to the minimum energy. Figure 3 shows the 
variation of the energy in two seperate cases. The first case utihses only the 
damping given by (15), the second curve uses a combination of the two. The 
differences are quite astounding. 
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Figure 3: Variation of the system energy due to damping. The first case 
involves a simple damping of the equations of motion as in (15) with v = 1.0. 
The second case involves a mixture of damping from the exciton field, and 
damping of the kinetic energy with v = 0.4. The magnitude of the exciton 
damping was set proportional to the derivative of the energy with respect to 
time. 
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Results showed that the quickest path to the stationary configuration 
involved making the exciton damping coefficient proportional to the gradient 
of the energy with respect to time. It was found that the exciton damping 
becomes less effective when the energy approaches that of the stationary 
configuration. Under these circumstances the damping of the kinetic energy 
becomes more useful, particularly for smaller values of v. Throughout the 
remainder of the investigation, a value of v — 0.4 was used in damping the 
system kinetic energy. 

3.2 Minimum energy configurations as a function of c 
and 1 

3.2.1 Parameters and Method 

The aim was to determine the minimum energy configurations for various 
values of the couplings I and c. The couplings j and k were set equal to 
1. Initially the couplings I and c were set to zero, and the minimum en- 
ergy configuration was found using the methods of the previous section. In 
fact, for these values of the couplings, the minimum energy configuration 
corresponded to the completely delocalized state, i.e. where the probability 
density for the exciton is approximately the same at all lattice sites. The 
value of I was increased incremently, and the static configuaration was again 
found. After increasing I upto 14 the coupling strength was incremently re- 
duced back to zero, following a similar procedure. Upon reaching zero the 
value of c was increased slightly, and the process repeated upto a maximum 
value of c = 22. 

3.2.2 Results 

Figure 4 shows the maximum in the probability density iV'ijImaa; 
minimum energy configuration as a function of the couplings I and c. The 
energy corresponding to these configurations is displayed in figure 5. 

As already mentioned, the ground state at (c, I) = (0, 0) corresponded 
to the completely delocalized state. Figure 4 shows how the region where 
this state corresponds to the ground state extends out to cover a significant 
proportion of the phase diagram for values of c approximately lying between 
the limits and 17. 

The various field configurations that were observed are shown in figure 6. 
Keeping / set equal to zero and increasing c had no effect on the configuration 
of the ground state until c ~ 8.5, when the state changed to that shown in 
figure 6a. Here the probability density for the exciton is localised mainly 
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Figure 4: Maximum in the probability density for the minimum energy con- 
figuration as a function of the couplings I and c. 

within a narrow band that wraps around the smallest lattice dimension. 
Upon increasing c further it was found that this state is stable until c ~ 18. 
For values of c higher than this, the configuration of least energy corresponds 
to the polaronic type state (fig 6b), localised mainly around a single lattice 
site. As expected the precise degree of localisation of the polaronic state was 
found to vary, depending upon the values of the couplings. The maximum 
in the probability density varied from 0.713 (c=21, 1=3) < [i/^l^aa, < 0.998 
(c=17, 1=0.25). 

Increasing / from the completely delocalized state at (c, /) = (0, 0) resulted 
in an observation of a different set of states. The delocalized state is the 
ground state until 1^3. Above this threshold value, the ground state once 
again corresponded to the wrapped state, already discussed. Increasing / 
still further resulted in the observation of states locahsed around six (fig 6c) 
and then four (fig 6d) sites respectively. Upon increasing c from this area 
of the phase diagram, the system underwent a smooth transition from that 
state mainly localised around four sites to a state principally surrounding 
two lattice sites. Taking / = 13.5 as an example, as the value of c increases 



12 



15 J. 

Figure 5: Energy of the ground state configuration as a function of the cou- 
plings I and c. 

over the region from 1 — 1.5, the probabihty density changed smoothly from 
that shown in fig 6d to that in fig 6f. The two-site state is then the ground 
state for c < 22, as is expected to be so for greater values of c. 

Two other configurations were observed (figure 6 g,h), however the states 
were always unstable. The state shown in fig 6g shows a two peaked state 
that is periodic around the nanotubc circumference. Fig 6h shows a kind of 
variation on the singly wrapped state of fig 6a. In this case there exist two 
separated bands of probability density wrapping the nanotube circumference. 

Each stationary state produced a distinct pattern of lattice displacement. 
In each case, the lattice sites are drawn in towards the points of high proba- 
bility density for the exciton. The characteristic pattern of distortion for the 
polaronic like state is depicted in figure 7a. The displacement of the lattice 
sites was considerable for those points located close to the centre of the prob- 
ability distribution. This could have important applications with respect to 
the electrical properties of hexagonal lattices, and therefore in carbon nan- 
otubes. An extremely interesting case is the wrapped state (fig 6a). The 
pattern of lattice distortion for this state is given by fig 7b. The transverse 
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displacement (perpendicular to the nanotube axis) is extremely small, ap- 
proximately of the order 10"^ the average longitudinal displacement. This 
state could therefore serve as a starting point for investigation into the effect 

of excitons upon the electrical properties of carbon nanotubcs. The state 
wraps around the circumference of the nanotube, and so any charge flowing 
parallel to the nanotube axis encounters approximately the same distortion. 

3.3 Stability Analysis of States 

The states discussed in the previous section were found to exist over larger 
areas of the phase diagram than Figure 4 suggests. Of course, in those regions 
they are higher energy unstable or metastable states. Thus when we perturb 
them sufficiently the system should fall to the ground state. This section 
details the results of perturbing states for values of c and I where they are 
not the minimum energy configuration. 

3.3.1 Polaronic States 

The method of perturbing the single-site localised polaronic-like states in- 
volved increasing the probability density marginally on the three nearest 
neighbour lattice sites and renormalising the exciton field. In order to quan- 
tify the stability of the states, perturbations of increasing strength were ap- 
plied and the strength of perturbation required to produce a transition was 
recorded. The perturbation parameter p specifies the strength. A strength 
of p = indicates no perturbation, while p = 1 corresponds to the three 
nearest neighbour sites being given equal probability density for the exciton 
as the initial central site. No perturbations with a strength greater than 1 
were induced. Figure 9 shows the strength of perturbation required to induce 
a transition as a function of the couplings / and c. For each distinct value of 
the couplings the initial conditions were given as those corresponding to the 
closest stable polaronic state. The value of p was then varied from to 1 in 
increments of 0.02. A strength of p = 1 indicates that the polaron state was 
stable for all applied perturbations, while a strength of p = means that the 
state decayed away naturally without any perturbation. For values of p in 
between, some degree of perturbation was required to knock the polaron into 
a lower energy state. This suggests that the state is a local, but not global, 
minimum in the corresponding region of the coupling phase space. 

Most values of the couplings displayed by the figure gave either stable 
(up to p = 1) polarons or completely unstable polarons. There are are some 
small regions where metastable states existed. These lie predominantly on 
the stability boundary between the polaronic states and the two-site localised 
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states. Metastable polaronic states were found to exist for values of / up to 
and including 1 — 5.25. For values of I greater than this the polaron state was 
found to naturally decay away to the two-site locahsed state. There exists 
another region of meta-stability between 3 < / < 5. This extends to smaller 
values of c < 17. In this region polaron states decayed to the completely 
delocalised configuration when subject to sufficient perturbation. 

3.3.2 Two-site-localised States 

A similar procedure was undertaken in order to perturb the states localised 
mainly around two lattice sites (figure 6f). The stability region of these states 
covers a large proportion of the phase diagram for c and I (figure 4). The 
initial conditions were those corresponding to the closest two-site localised 
state. For each value of c the value of / was decreased by different amounts 
and perturbations were applied, again the aim being to cause a spreading 
out of the sohton into one of the lower energy states. The perturbations 
involved increasing the exciton field slightly at the four nearest neighbour 
sites and renormalising the field. The level of perturbation required to induce 
a transition was investigated. 

The perturbations were in this case ineffective in causing transitions. The 
states were stable up to perturbations of magnitude p — 1. However for low 
enough values of I the state decayed naturally into a lower energy state: 
either a singly-wrapped state (figure 6a) or the completely delocalised state. 
Figure 10 shows the maximum in the probability density as / was decreased 
for different values of c. The height of the central peak on the two sites was 
found to decrease with /. To the left of the black fine the minimum energy 
state is that with the exciton field localised mainly around two sites. To 
the right of this line the minimum energy state corresponded either to the 
singly- wrapped state or the entirely delocalised state (see figure 4). It is clear 
from the figure that there exists a region where the states localised around 
two lattice sites show a large degree of stability, even though the minimum 
energy configuration was some other state. Furthermore, this region covers 
a significant proportion of the / versus c phase diagram. 

3.3.3 Stability of other states 

Perturbations were applied to some of the other states displayed in figure 
6. The singly- wrapped state (figure 6a) was studied in this way: in particu- 
lar the boundary between the minimum energy phases of the singly-wrapped 
state and the completely delocalised state starting at (c, /) ^ (8.75, 0). Start- 
ing at these values of the couplings, I was gradually increased. The singly- 
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wrapped state was found to decay into the completely delocalised state with 
extremely little or indeed no perturbation. In each case the transition was 
extremely gradual, with the maximum in the probability density slowly de- 
creasing with time. This was in stark contrast to the transitions observed 
from the polaronic and two-site localised states, which were found to decay 
extremely quickly. The procedure was repeated along the boundary, giving 
similar results. The state was found to show a large degree of stability only 
within those regions where it was the minimum energy state. However this 
is a considerable part of the c, I phase diagram. 

The stability of the state localised around six sites was also studied in 
detail for decreasing /. In this case perturbations were applied once again 
by increasing the magnitude of the exciton field around the soliton in an 
attempt to encourage it to spread out. Figure 11 displays the results. 

The region to the left of the circles corresponds to the region in figure 
4 where the minimum energy state is the six-site-localised state. To the 
right of the circles the minimum energy configuration is the singly-wrapped 
state. Perturbations were again applied by increasing the exciton field at the 
nearest neighbour sites. The narrow band between the hues of circles and 
crosses in figure 11 is the region where the six-site solitons showed a degree of 
stability. In fact the states were stable up to values of p = 1 for all values of 
the couplings / and c in this region. In the region to the right of the crosses 
the solitons localised around six sites were found to be unstable, decaying 
into the singly-wrapped state in the absence of any applied perturbation. 

4 Conclusion 

Terms describing the coupling of the exciton resonance interaction to lattice 
displacements were included into a Hamiltonian describing the existence of 
excitonic states on a hexagonal lattice with continuous boundary conditions. 
The added terms were found to spoil the elegant analj^ic solution to the equa- 
tions of motion found in [1] , however numerical simulation showed that these 
terms produced a rich structure in the coupling phase space of stationary 
states. By damping the displacement and exciton fields, solitonic stationary 
states localised around varying numbers of lattice sites were observed includ- 
ing, in agreement with previous results, extremely localised polaron states. 
Certain solitons were also found to wrap around the lattice, producing pat- 
terns of lattice displacement that may perhaps be used to study the effects 
of solitons upon the electrical properties of carbon nanotubes. The stability 
of states for certain values of c and / was investigated numerically by intro- 
ducing perturbations. In most cases the system possessed metastable states 
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in addition to the ground state. Some of these metastable states appeared 
to be separated from the ground state by a small energy barrier and were 
thus local minima of the energy. Other excited states were found to exhibit 
stability under perturbation. 



5 References 

[1] B. Hartmann and W.J. Zakrzewski, Electrons on Hexagonal Lattices and 

applications to nanotubes, Phys.Rev. B 68 (2003) 184302 

[2] A.S.Davydov and N.I. Kislukha 1973, Phys. Stat. Sol. (b) 59,456 

[3] A.S. Davydov, Solitons in Molecular Systems 2nd e(i.,Kluwer Academic 

(1991) 

[4] L.Brizhik, A. Eremko, B. Piette and W.J. Zakrzewski, Physica D 146 

(2000) , 275. 

[5] L.Brizhik, B. Piette and W.J. Zakrzewski, Ukr. Fiz. Journal 46 (2001), 
503. 

[6] L.Brizhik, A. Eremko, B. Piette and W.J. Zakrzewski, Physica D 159 

(2001) , 71. 

[7] J.C. Eilbeck and M. Johansson, The Discrete Nonlinear Scroedinger Equa- 
tion - 20 years on, |arXiv:nlin/0211049| 



[8] P.G. Kevrekidis, B.A. Malomed and Yu. B. Gaididei, Solitons in triangu- 
lar and honeycomb dynamical lattices with the cubic nonlinearity, Phys. Rev. 
E 66 (2002) 016609 

[9] R. Saito, G. Dresselhaus and M. S. Dresselhaus, Physical properties of 
carbon nanotubes, World Scientific (1998) 



6 Appendix Al - Hamiltonian and equations 
of motion for simulation 

The simulated equations are altered to account for four different lattice sites 
separately. The advantage of this method lies in the reduction of the array 
sizes used throughout the simulations. For more information see [1]. The 
Hamiltonian used throughout was: 



2 2 



^ = E E iE + W)^P,,^Pl^-J^Pl^i^P,+,J+^^.,,+^P,+,J^^^ (20) 
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The corresponding equations of motion are, for ip: 
i = l + Ak 
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while equations describing the lattice site displacements Ui^j and Vjj are: 
i = 1 + 4fc 
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-I 



i^2 + Ak 



i^3 + Ak 
(Pu 



k{3uij — Ui+ij — Mj-ij — Mj-ij+i) 



3M 



(25) 



(26) 
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^ — k{3uij — itj+ij — Ui-ij — 



and i = 1 + Ak 



2 + Ak 



3 + Ak 



d^Vi 



df^ 



d'^Vi 



dt^ 



— (2|^i+i,jf - l-^i-ij+il^ - iV'i-ijI 



(l^i+ij'T - IV^i+ij-iT) 



-iV3 



+/\/3 



\/3M 



+lV3 
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k{3vi 



Vi-l,j - Vi_ij 
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Figure 6: Probability density distribution over the hexagonal lattice of the 
observed states. 




Figure 7: Pattern of lattice displacement centred around a polaronic type 
state. 
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;ure 8: Pattern of lattice displacement for the singly wrapped state. 
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Figure 9: Perturbation strength required to induce a transition from the po- 
laronic state as a function of c and I 
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Figure 10: Stability of the two-site-localised state as a function of c and I 
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Figure 11: Stability of the six- site- localised state as a function of c and I 
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